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This  paper  develops  an  algorithm  for  recovering  a  collection  of  linear  cracks  in  a  homogeneous 
electrical  conductor  from  boundary  measurements  of  voltages  induced  by  specified  current 
fluxes.  The  technique  is  a  variation  of  Newton’s  method  and  is  based  on  taking  weighted 
averages  of  the  boundary  data.  The  method  also  adaptively  changes  the  applied  current  flux 
at  each  iteration  to  maintain  maximum  sensitivity  to  the  estimated  locations  of  the  cracks. 
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1  Introduction 


In  this  paper  we  develop  a  very  efficient  computational  algorithm  to  reconstruct  a  collection 
of  linear  cracks  inside  a  homogeneous  conductor  from  electrostatic  boundary  measurements. 
The  algorithm  in  this  paper  can  be  seen  as  a  natural  extension  of  the  algorithm  developed  in 
[15]  for  the  reconstruction  of  a  single  crack.  This  extension  poses  several  theoretical  as  well 
as  practical  challenges.  We  have  also  significantly  improved  the  efficiency  and  versatility  of 
the  earlier  algorithm  by  basing  all  computations  for  the  underlying  conductance  problem  on 
a  one-dimensional  boundary  integral  formulation  instead  of  a  two-dimensional  finite  element 
formulation.  It  should  be  mentioned  here  that  boundary  integral  formulations  have  been  used 
in  other  implementations  of  (single)  crack  reconstruction  algorithms  [13],  [14],  and  also  that 
progress  on  the  development  of  an  algorithm  for  the  reconstruction  of  a  single  (penny-shaped) 
crack  inside  a  three  dimensional  object  is  reported  in  [14].  Algorithms  like  the  one  discussed 
in  this  paper  are  significantly  different  from  more  general  purpose  imaging  algorithms  (cf. 
[3],  [6],  [9],  [16])  that  seek  to  reconstruct  an  unknown  distributed  conductivity  profile  from 
similar  boundary  measurements.  Algorithms  such  as  ours  are  based  on  the  assumption  that 
certain  apriori  information  about  the  profile  is  available  and  they  incorporate  this  knowledge 
into  the  reconstruction  in  such  a  way  as  to  achieve  better  continuous  dependence  and  better 
discrete  approximation  properties.  One  important  feature  of  the  present  algorithm  is  that 
it  is  based  on  an  adaptive  change  of  the  prescribed  boundary  current  patterns  to  ensure 
“maximal”  sensitivity.  The  idea  to  use  some  kind  of  “optimal”  current  pattern  in  connection 
with  impedance  imaging  has  been  developed  by  Gisser,  Isaacson  and  Newell  (cf.  [9]);  the 
specific  strategy  we  use  is  somewhat  different  from  theirs  and  ties  in  directly  with  the  iterative 
procedure  (it  does  not  rely  on  any  eigenfunctions).  Our  reconstruction  is  based  on  the  usage 
of  relatively  few  averages  of  the  boundary  voltage  measurements  (as  opposed  to  all  the 
boundary  voltage  data).  In  addition  to  improving  the  efficiency  of  the  algorithm,  this  should 
also  decrease  the  probability  of  getting  caught  in  a  local  minimum  when  compared  to  more 
standard  output  least-squares  algorithms. 

An  outline  of  this  paper  is  as  follows.  In  Section  2  we  present  the  “customary”  mathe¬ 
matical  model  of  electrostatic  conductance  for  a  smooth,  isotropic  background  medium  that 
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contains  a  collection  of  cracks.  In  particular  we  demonstrate  the  duality  between  the  notions 
of  perfectly  insulating  cracks  and  perfectly  conducting  cracks.  We  also  briefly  discuss  known 
uniqueness  and  continuous  dependence  results.  Section  3  contains  a  detailed  description  of 
the  4 n  functionals  that  we  use  for  the  reconstruction  of  a  collection  of  n  or  fewer  linear  cracks. 
This  section  also  provides  a  discussion  of  the  adaptive  strategy  that  we  use  for  the  selection 
of  the  “maximally  sensitive”  electrode  locations.  The  boundary  integral  formulation  of  the 
electrostatic  conductance  problem  and  its  discretization  by  Nystrom’s  method  is  the  topic  of 
Section  4.  The  centra]  part  of  the  reconstruction  algorithm  is  a  version  of  Newton’s  method. 
This  particular  version  together  with  the  required  gradient  computation  is  discussed  in  Sec¬ 
tion  5.  Section  6  contains  a  selection  of  representative  computational  experiments  with  our 
algorithm.  Finally  in  Section  7  we  provide  a  brief  summary  of  our  results  and  a  description 
of  possible  future  developments. 


2  The  Mathematical  Model 

A  single  crack  is  commonly  modeled  as  a  perfectly  insulating  curve  a.  With  a  background 
conductivity  0  <  70  <  7(2)  <  71  and  a  finite  collection  of  cracks  E  =  U£=1<Tjt,  the  steady 
state  conductance  equations  thus  read 

V  •  (jVv)  =  0  in  0  \  E,  (2.1) 

dv 

1dv  =  0  on 

with  appropriate  boundary  conditions  on  dfl,  e.g., 

v  =  <f>  on  dfl.  (2.2) 

The  field  v  is  normal  to  E.  The  function  v  represents  the  voltage  potential  induced  in  0. 
We  assume  that  0  is  simply  connected,  i.e.,  it  has  no  holes,  and  so  the  entire  boundary  dfl 
is  accessible  from  the  “outside”.  Let  u  denote  the  “7-harmonic”  conjugate  to  v.  It  is  related 
to  v  by  the  formula 

(Vu)x  =  yVv,  (2.3) 
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where  X  indicates  counter-clockwise  rotation  by  ir/2.  For  a  particular  set  of  constants 
Cfc,  k  =  1,  ...  ,n,  the  function  u  solves  the  problem 


V  •  (7”1  Vu)  =  0  in  ft\E, 


u  =  Ck  on  <7fe,  fc  =  1,  ...  ,n 


d<l>  „ 

7  ar  =  *  =  aJ on  aa 


Here  s  denotes  the  counter-clockwise  tangent  direction  on  5ft  and  v  denotes  the  outward 
normal  on  5ft.  For  these  particular  constants,  finding  a  solution  to  (2.1),  (2.2)  is  thus 
equivalent  to  finding  a  solution  to  (2.4),  (2.5).  The  constants  c*  may  (up  to  a  common 
additive  constant)  be  characterized  in  several  equivalent  ways: 


(a)  Let  A  be  the  n  x  n-matrix  with  elements  Aij  =  fa  y~iVU^VU^J^dx  and  let  b  be  the 
n-vector  with  elements  6,-  =  fan  t/iU^ds,  where  U^\  j  =  1, . . . ,  n,  denote  the  solutions 
to 


V  (7 =  0  in  ft\E, 

=  1  on  <7 j, 

=  0  on  cr/e,  k  ^  j 


with 


Then  the  vector  c  =  (ci, 


(b)  The  ^-harmonic”  conjugate,  u,  satisfies  fVk[ 7-1f^]  ds  =  0,  1  <  fc  <  n.  Here 

[7~ 1  a^3  =  7_1^r —  7-1^  denotes  the  jump  in  the  normal  flux  across  the  curve  <r*.  2 
Furthermore,  the  set  of  constants  {c*}^  is  the  unique  set  of  constants  for  which  the 
solution  to  (2.4),  (2.5)  has  this  property. 

3The  expression  -§~  denotes  the  limit  of  the  derivative  (in  the  direction  v)  as  one  approaches  <r*  from 
the  side  to  which  v  points.  denotes  the  limit  as  one  approaches  <r*  from  the  opposite  side. 
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(c)  Let  T  be  a  fixed  point  on  dfi,  in  a  neighborhood  of  which  <f>  is  smooth.  Let  r*  be 
a  smooth  curve  in  il  \  S  connecting  T  to  an  interior  point  of  the  crack  and  let 
s  denote  the  unit  tangent  direction  along  r*,  pointing  from  T  towards  <7*.  Then  the 
constants  c*  are  given  by  the  formulae 

c‘  =  -/.7  %ia+u^ 
where  v  denotes  the  normal  field  u  =  — sx. 

The  characterization  (c)  is  a  direct  consequence  of  the  relation  (2.3).  The  characteri¬ 
zations  (a)  and  (b)  are  practically  much  more  useful;  they  are  both  a  consequence  of  the 
following  well  known  result  from  convex  duality. 

Proposition  2.1  If  <f>  is  an  element  of  Hx^2(dCl),  then  the  field  £  =  q/'Vv  is  the  (unique) 
minimizer  of  the  functional 

1  dx  ~  Ln^V’t/ds  (2-6) 

2  Jn  Jan 

in  the  space  H  =  L2(Cl)  D  {q  :  V  •  t]  =  0  in  Cl  \  £  ,  rj  •  v  =  0  on  <7*,  k=l,. . .  ,n}. 

It  is  not  difficult  to  see  that  any  element  of  the  space  H  satisfies  V  •  q  —  0  in  all  of 
ft  and  therefore  has  the  form  q  =  (Vu>)x  for  some  w  6  /^(ft),  with  w  being  constant  on 
each  <7*.  Conversely,  it  is  also  true  that  any  vector  field  of  the  form  rj  =  (Vtn)x,  w  €  //'(ft) 
H{u)  =  constant  on  each  <7 k  =  1,  ...  ,n},  is  in  the  space  H.  After  insertion  into  (2.6)  we 
obtain  7VU  =  (Vu)x,  where  u  is  the  minimizer  of  the  functional 

^  /  7_1|Vu>|2  dx  -  I  —w  ds 
2Jn  Jan  os 

in  the  space  //'(ft) H  {w  =  constant  on  each  <r*,  k  =  1,  ...  ,n}.  This  provides  a  variational 
characterization  of  the  ^-harmonic”  conjugate  to  v.  Let  F(d),  d  €  IRn,  denote  the  expression 

m  =  U  (  «..+£>(/<*>)<(*• 

2  Jn  fc=l  Jdn  k^l 

From  the  above  discussion  it  is  clear  that  the  set  of  constants  corresponding  to  u  are  char¬ 
acterized  by  the  fact  that  d  =  0  is  a  minimum  of  F.  Equivalently  (because  of  the  form  of 
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F),  the  set  of  constants  corresponding  to  u  are  characterized  by  the  fact  that  d  =  0  is  a 
stationary  point  for  F.  Stationarity  of  d  =  0  is  equivalent  to  the  conditions 


/  7-1VuVf/(*)  dx=  I  xl>U(k)  ds,  k  =  1,  ...  ,  n.  (2.7) 

J  n  ./an 

From  integration  by  parts  (on  the  domain  fi  \  E)  we  have 

/  7-1VuVf/(fc)  dx=  I  [7-1^]  ds  +  /  <k. 

7n  Ja  n 

Insertion  of  this  into  (2.7)  now  gives  that  the  set  of  constants  corresponding  to  u  are  char¬ 
acterized  by 

L  =  0  k  =  \,  ...  , n, 

as  asserted  in  (b). 

On  the  other  hand,  the  function  u  has  the  form 

u  =  uo  +  (2.8) 

«=i 


where  u0  is  the  solution  to 


V-(7-xVu0)  =  0  in  fi\E, 


uq  =  0  on  <7j t,  fc  =  l,...,n, 


with 

7-1^^  =  xp  on  50.  (2-10) 

au 

Because  of  (2.9)  and  (2.10)  we  have  /n  7_1Vu0V{/^l  dx  =  0;  by  insertion  of  (2.8)  into  (2.7) 
and  use  of  this  formula  we  now  obtain 

fcj/  7-1V(/(,)W<*>  dx  =  f  xl>U{k)  ds  k  =  l,...,n, 

/o  ./an 

which  is  exactly  the  characterization  (a).  The  above  argument  rests  on  the  fact  that  <p  is 
in  H1^2(dil)  (xp  is  in  //-1/f2(50)).  However,  by  continuity  the  characterizations  carry  over 
to  cases  in  which  <f>  (and  xp)  are  not  necessarily  so  regular  that  u  is  a  variational  solution. 
In  particular  these  characterizations  remain  valid  if  xp  consists  of  delta  functions,  a  type  of 
boundary  current  we  shall  repeatedly  use  in  this  paper. 


5 


Since  u  and  v  are  related  by  the  equation  ( Vu )J-  =  7V1;,  it  is  clear  that  knowledge  of  the 
pair  (<f>,  7§^|an)  is  equivalent  to  knowledge  of  the  pair  (u|an,  ip).  It  is  much  more  convenient 
to  work  with  the  function  u  as  opposed  to  the  function  v,  and  we  shall  entirely  do  so  for 
the  development  of  our  algorithm.  In  particular  by  working  with  u  we  avoid  the  difficulties 
that  are  associated  with  non  integrable  kernels  in  the  integral  equation  formulation  (cf.  [12], 

[14])- 

Let  Pq ,  ...  ,  Pm  be  M  +  1  points  on  dfl;  we  assume  that  these  points  are  labeled  in  order 
of  counter-clockwise  appearance,  starting  from  Pq.  For  the  crack  reconstruction  we  utilize 


solutions  corresponding  to  the  two-electrode  currents  ipj  =  8Po  -8P},  j  =  1,...,M, 


V  •  (7  1  Vuj)  =  0  in  fl  \  E, 

Uj  =  cjjfl  on  (Tk,  k  =  1,  ...  ,n, 


(2.11) 


with 

/)•»/  • 

=  <p,-«p,  on  an,  (2.12) 

the  constants  cj^  being  selected  so  that  Xrfc[7-1^]  ds  =  0  for  &  =  1,  ...  ,n.  The  inverse 
problem  may  now  be  stated  explicitly  as  follows: 

We  seek  to  reconstruct  the  collection  of  cracks  E  =  UJLjO*  from  knowledge  of  the  boundary 
voltage  data  {ujlanjjii  corresponding  to  the  prescribed  two-electrode  currents  7_1  ^  =  8p0  —  6p- , 
j  =  1,  ...  ,  M . 

It  is  known  that  boundary  voltage  measurements  corresponding  to  M  =  n  +  1  fixed  two- 
electrode  currents  suffice  to  uniquely  identify  a  collection  of  n  (or  fewer)  cracks  [4].  This 
result  is  an  extension  of  a  result  in  [8]  which  asserts  that  boundary  voltage  measurements 
corresponding  to  two  fixed  two-electrode  currents  suffice  to  uniquely  identify  a  single  crack. 
Recently  an  interesting  continuous  dependence  estimate  has  been  obtained  for  the  case  when 
the  background  conductivity  is  constant  and  there  is  at  most  one  crack  [1].  Briefly  described, 
this  estimate  states  that  if  the  boundary  voltage  data  (on  some  open  subset  of  dil)  deviate 
by  e  then  the  crack  locations  differ  by  at  most  [log(|  log  e|)]-1/4.  In  the  present  paper  shall 
always  try  to  fit  the  data  entireiy  by  means  of  linear  cracks.  For  such  cracks  one  can 
hope  to  have  better  continuous  dependence  estimates,  as  indicated  by  the  results  in  [8].  As 
was  the  case  in  [15],  we  base  the  reconstruction  of  the  cracks  on  the  values  of  a  number 
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of  functionals  (as  opposed  to  all  the  boundary  voltage  measurements).  In  [15]  we  used  4 
functionals  corresponding  to  the  reconstruction  of  a  single  linear  crack;  the  natural  extension 
is  to  use  4n  functionals  for  the  reconstruction  of  n  cracks.  In  the  following  section  we  give 
a  careful  description  of  these  functionals. 


3  The  Functionals 

We  now  specialize  to  the  case  of  a  constant  background  conductivity,  7=1.  Let  F  denote 
the  vector-valued  function 

F(£,  A  w)  -  (F(£,  »">),  F(£,  F( E,  V>,  u,<3>),  F(E,  <A,  uj'4’))', 

where  F(T,,ip,w)  is  given  by 

F(Z,i>,w)=  [  (3.1) 

and  where  1  <  i  <  4,  are  particular  solutions  of 

Au;  =  0  in  IR2  \  E. 


The  function  u  =  is  the  solution  to 


A u  =  0  in  n\S, 


u 

du 

dv 


ck  on  (Tk ,  A;  =  l,...,n, 
ip  on  dfi, 


with  the  constants  c*  uniquely  specified  by 


uds  =  0 


and 


ds  =  0, 


k  =  1, . . 


n. 


(3.2) 


The  exact  selection  of  boundary  currents  ip  and  test  functions  w  =  (w^\w^2\w^3\w^Y  is 
very  important  and  will  be  discussed  shortly.  We  select  one  ip  and  one  w  corresponding  to 
each  crack  ak',  whenever  we  want  to  emphasize  this  correspondence  we  use  the  notation 


rp 


<?k 


and 
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We  have  for  convenience  chosen  the  normalization  fdn  u  ds  =  0  for  the  voltage  potential.  We 
shall  always  select  w  so  that 

f  dw 1*1  r  rau»Wl 

/  —5 — ds  =  0,  and  /  — —  ds  =  0  Ar  =  l,...,n.  (3.3) 

Ja n  ai/  dv 

Because  of  the  first  identity  in  (3.3)  the  function  F  is  unchanged  by  the  addition  of  a 
constant  to  u,  and  we  could  therefore  just  as  easily  work  with  any  other  normalization. 
The  components  of  F  are  just  weighted  averages  of  the  boundary  voltage  data.  We  use  a 
weighting  function  of  the  form  because  of  the  relation 

[  dw  t  * 

/  tt(E,0)-r—  ds  =  /  Vii(E,0)Vtc  dx  (3.4) 

Jan  d  v  Jn\z 

that  exists  between  the  expression  in  equation  (3.1)  and  the  energy  bilinear  form  (by  means 
of  Green’s  formula).  As  will  be  seen  later,  these  averages  are  equivalent,  in  the  absence  of 
any  crack,  to  the  set  of  the  first  4  nontrivial  Fourier  modes  of  the  induced  boundary  voltage. 

The  data  for  our  reconstruction  consist  of  measured  boundary  potentials  corresponding 
to  certain  prescribed  two-electrode  boundary  currents.  We  denote  the  voltage  data  cor¬ 
responding  to  the  boundary  current  ip  by  g(ip)  and  define  a  corresponding  vector-valued 
function 

f  (0,  w)  =  (f(rp,w{1)),f(ip,w{2)),f(ip,w(3))J(ip,w(4))y, 
where  f(ip,  w )  is  given  by 

f(ip,w)  =  f  g(tp)~  ds, 

Jan  ov 

and  where  itn*l  are  the  same  functions  as  before.  Our  algorithm  seeks  a  solution  S  = 
to  the  4n  equations 

F(S,0<r*,wffik)  =  f  (0<r*,wffJ,  1  <  fc  <  n. 

Consequently,  we  do  not  use  information  about  the  full  boundary  voltages  for  the  recon¬ 
struction;  we  only  use  information  about  the  values  of  these  particular  functionals. 

We  implicitly  assume  that  our  data  is  consistent  so  that  f  {ip„ki  w^)  corresponds  tu  some 
collection  of  cracks  E*  (E*  may  consist  of  fewer  than  n  cracks).  In  reality  we  therefore  solve 

G*(E)  =  0,  \<k<n ,  (3.5) 
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with 


G*(E)  =  F(E,^„w„)-F(S-,^,w„,)  (3.6) 

Clearly  E*  is  a  solution  to  (3.5).  If  the  Frechet  derivative  /?s{Gjt}|s=s*  (a  4nx4n  matrix) 
is  nonsingular,  then  E*  is  indeed  the  unique  solution  to  (3.5)  near  E*,  and  furthermore,  one 
may  expect  that  some  variation  of  Newton’s  method  will  be  an  efficient  solution  technique. 
Differentiation  with  application  of  the  “chain  rule”  yields  the  expression 

D£{Gfc}£=1|s=s.  =  {DsF(E,V’<Tj,w(Tj)|E=s.}fc=1, 

for  the  derivative  with  respect  to  E  (at  E*),  the  right  hand  side  of  which  is  a  4n  x  4n  matrix 
with  rows 

i \  9u'L* 

DeF(E,0<t.,uj^)|s=£.  =  fanD su(E,V'<Tj;)|E=E*-a^L  ds,  (3.7) 

1  <  k  <  n,  1  <  *  <  4. 

In  [15]  we  explicitly  calculated  the  expression  (3.7)  in  terms  of  u  and  (eliminating  D^u); 
we  used  this  alternate  expression  for  the  selection  of  the  “maximally  sensitive  two-electrode” 
currents  as  well  as  for  the  Newton’s  update  itself.  In  our  implementation  here  we  shall  rely 
on  essentially  the  same  technique  as  in  [15]  for  the  selection  of  two-electrode  currents,  but 
we  shall  directly  compute  the  derivative  of  u(E,t/>)  with  respect  to  E  at  the  discrete  level 
(c/.  Section  5)  for  use  with  the  Newton’s  update. 

When  talking  about  the  derivative  with  respect  to  E,  we  mean  the  derivative  with  respect 
to  the  4n  parameters  that  are  used  to  describe  E.  Just  as  in  [15]  we  parametrize  a  single 
crack  by  (b i,  63),  0,  and  A,  where  (f>i,  62)  are  the  coordinates  of  one  endpoint,  9  is  the  counter¬ 
clockwise  angle  between  the  crack  and  the  halfline  y  =  &2,  x  >  b\  ,  and  A  is  the  length  of  the 
crack  (here  coordinates  of  points  are  relative  to  a  fixed  reference  coordinate  system).  For 
convenience  we  order  these  coordinates  q  =  ( 9 ,  62, 61,  A);  the  parameter  set  corresponding  to 
E  is  now  given  by  the  vector  Q  =  (qi,q2,  ...  ,  qn)*-  The  derivative  Dz{Gk}  consists  of  the 
(n3)  4  x  4  blocks 

DqiGk. 
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At  the  solution,  £  =  £*,  these  blocks  equal 


0qlF(Q,ifc,.,w4,.)|Q=Q..  (3.8) 

If  the  derivatives  (3.8)  are  formed  at  some  Q°  not  equal  to  Q*  then  they  no  longer  represent 
the  full  derivative  of  {G*};  the  latter  also  includes  terms  where  and  \\)„k  are  differentiated 
through  (c/.  (3.6)).  For  the  Newton’s  update  we  use  a  more  complete  “derivative”  which 
includes  the  terms  that  arise  when  the  wtTk  are  differentiated  through,  but  we  do  not  include 
differentiation  through  Vv*-  The  goal  behind  the  choice  of  w,k  and  xpak  is  to  make  the 
derivative  of  {Gfc  }  as  far  from  singular  as  possible  at  E  =?  E°,  the  current  stage  of  the 
iteration. 

In  order  to  describe  the  choice  of  we  select  a  coordinate  system  such  that  er^  lies  on 
the  positive  X\  axis,  with  one  endpoint  at  the  origin.  In  this  coordinate  system  we  choose 


wi',’  =  Im(z],  to™  =  Im[z2], 


(3.9) 


<31  l  Ml2  -  **)/*(* -**)].  Re(z)  >  ¥  , 

<’  =  {  (3.10) 

(  -Re[(z  -  A,) y/z(z  -  A,)],  Re(z)  < 


Re(\/z(z  -  •'*))' 


Re(z)  >  ^ 


(3.11) 


[  -R e[yjz(z  -  A*)],  Re(z)  <  **■ 

where  z  =  x j  -\-ix2  and  A*  denotes  the  length  of  <7*.  The  functions  and  are  extended 
to  Re(z)  =  Afc/2  by  continuity.  Except  for  a  change  in  this  is  exactly  the  same  choice 
of  test  functions  as  in  [15].  Since  they  are  harmonic  in  fi,  the  two  functions  wy^  and 
clearly  satisfy  (3.3).  It  requires  some  extra  calculations  to  check  that  and  wyj  also 
satisfy  (3.3);  for  reasons  of  brevity  we  omit  these  calculations  here.  The  fact  that  w^3J  and 
do  indeed  satisfy  (3.3)  makes  the  Remark  1  in  Section  3  of  [15]  superfluous.  Notice  that 
w^,  1  5:  *  <  4,  vanish  on  the  crack  cr*. 

Remark 

Given  the  specific  form  of  the  weight  functions  it  is  now  fairly  easy  to  explain  why,  in 
the  definition  of  G*,  we  pick  a  distinct  boundary  current  t/v*  corresponding  to  each  k.  If  as 
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an  extreme  we  had  picked  the  same  boundary  current  corresponding  to  each  crack,  then 
the  first  two  equations  of  (3.5),  (3.6)  would  be  solved  for  all  k  if 


Lb(y)u(Tj'i,)  is 

Li;ix2-y,)u(^)ds 

f  d 

I  «-(*y)u(£,V0  ds 

Je  n  a  v 


=  Lrv(x)9Wia' 

= 

=  L  ry  - y2)gw  ds 

=  /  «-( xy)gW 

Jdil  OV 


where  ( x,y )  denote  coordinates  relative  to  some  fixed  coordinate  system.  The  system  (3.5), 
(3.6)  would  therefore  represent  1. 1  more  than  2n  +  4  equations  for  the  4 n  unknowns  of  E. 
For  n  >  2  this  immediately  leads  to  “underdetermination”  and  a  singular  Jacobian.  ■ 

In  [15]  we  analyzed  the  structure  of  the  derivative 


^qF(*l>  W(r°)|q=q°>  (3.12) 

for  the  case  of  a  single  crack  (and  with  a  slightly  different  choice  of  u^.4)).  With  the  ordering 
of  the  parameters  ( 6 ,  b2 ,  b\,  A)  we  found  that  this  4x4  matrix  was  lower  triangular.  We  also 
found  that  if  the  test  functions  it/*),  1  <  i  <  4,  had  been  selected  harmonic  (and  0  on  the 
crack)  then  the  last  two  columns  in  this  matrix  would  have  been  zero.  Test  functions  with 
singularities  like  it/3)  and  it/4)  are  thus  essential  to  insure  that  this  matrix  is  non-singular. 
Since  the  only  change  we  make  in  the  test  functions  concern  u/4)  we  do  not  destroy  the  lower 
triangular  structure  of  this  derivative  for  the  one  crack  case.  In  the  multiple  crack  case  the 
counterparts  of  the  matrix  just  discussed  are  the  diagonal  entries 


£>q*F(Q>^<705w<r°)lQ=Q0- 


It  is  worthwhile  noticing  that  these  matrices  do  not  inherit  the  lower  triangular  structure. 
For  a  any  fixed  crack  <7*,  contributions  corresponding  to  the  other  cracks  will  appear  above 
the  diagonal.  These  contributions  are  of  the  form  —  ds ,  where  the  functions  w^k  3 

are  related  to  w £’)  by  the  formulae 


w 


(’)  _ 


w 


=  (-X2,Xi  )‘*Vu>g,  117*2,2  = 


dx2 


(')  _ 
<*k  .3 


dw(f\ 

dx\  ' 


wl'1,4  =  (*I,*2)*- VwJJ, 


(3.13) 
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cf.  [15],  page  921.  However,  it  is  our  practical  experience  that  the  above  selection  of 
together  with  the  appropriate  selection  of  the  two-electrode  current  ipak  (to  be  discussed 
below)  is  a  very  effective  way  of  achieving  a  Frechet  derivative  which  is  far  from  singular. 

Remark 

It  is  interesting  to  consider  the  limit  A*  — ►  0  as  one  endpoint  and  the  direction  of  the  crack 
stay  fixed.  Let  crfc  be  a  crack  with  endpoint  at  (61,62),  zero  angle,  and  length  A/t-  The 
corresponding  limits  of  the  functions  w b)  are  given  in  absolute  coordinates  by 

wo}  =  y-b  2,  u;J2)  =  2(x  -  6i)(y  -  62),  (3.14) 

w, o3)  =  (x  -  6j)2  -  (y  -  62)2,  Wq4^  =  x  —  bj. 

Actually,  w^J  =  Wq ^  and  identically,  since  these  two  functions  do  not  depend 

on  the  crack  length.  Moreover,  the  approximations  Wq3^  and  «  w for  the  third 

and  fourth  functionals  are  quite  accurate  sufficiently  far  away  (e.g.  two  crack  lengths)  from 
the  crack.  If  fl  is  the  unit  ball  and  g(9),  0  <  9  <  2ir,  is  the  limiting  boundary  voltage,  and 
F°  =  limA^o  F(£,  V>,  tn^),  then  (3.14)  gives 

F®  =  /q  w  g{9)  sin  9  dO 

F$  =  2  Jo2"  g(0)  sin  29  d6  -  26,  g{9)  sin  9  dO  -  262  /02*  g(9)  cos  9  d9 

F3°  =  2  ft*  g(9)  cos  29  d9  -  26,  g(9)  cos  9  d9  +  262  g(9)  sin  9  d9 

F°  =  ft*  g(9)  cos  9  d9. 

In  the  limit  A*  — >  0,  F  thus  represents  the  first  4  Fourier  modes  of  the  boundary  data.  ■ 
The  change  we  have  made  in  can  be  explained  from  this  remark:  with  the  choice  made 
in  [15]  u/3)  and  had  the  same  limit  as  A  — ►  0.  The  boundary  integral  implementation  of 
the  algorithm  would  therefore  occasionally  attempt  to  fit  the  data  by  just  making  the  cracks 
very  small,  since  it  then  in  effect  only  would  have  to  satisfy  3n  equations  (using  3n  unknowns) 
as  opposed  to  satisfying  a  larger  set  of  4n  equations.  We  did  not  encounter  this  phenomenon 
in  any  of  the  single  crack  experiments  performed  in  [15]  since  the  implementation  there 
was  based  on  a  2-d  finite  element  formulation,  which  effectively  put  a  lower  bound  on  the 
crack  length.  The  new  t^4)  we  have  introduced  here  is  simply  a  (correctly  scaled)  linear 
combination  of  the  and  w ^  used  in  [15]. 
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For  the  selection  of  the  current  tj)0k  we  rely  on  an  appropriate  adaptation  of  the  technique 
developed  in  [15].  In  that  paper  we  calculated  that  the  first  and  second  diagonal  entries  of 
the  matrix  (3.12)  are  m  and  2m  respectively,  where  m  denotes  the  expression 

m  =  / 

Jan 

We  then  proceeded  to  select  a  two-electrode  current  which  made  this  expression  largest 
possible.  In  the  multiple  crack  case  the  first  diagonal  entry  of  the  k’th  diagonal  block  of 
{DEF(E,^fc,w£rJ}2=1  is  given  by 

mk = L  is  -  /,  [It]  ”«•< is  ’  <3-15> 

where  is  as  defined  by  (3.13). 

Consider  the  function  (k  €  Z/1  (fl  \  <?k)  satisfying 


ds, 

ov  1 


A  (k  =  0  in  Q,\crk, 


6  =  0  on  <xk ,  (3.16) 

dtk  dw^ 

—  =  ■■■  ■  on  Oil. 

Ov  ov 

Note  that  does  not  identically  vanish  on  crk,  so  that  £k  is  different  from  After 

insertion  of  £k  into  (3.15)  we  have 

/  (Mk  .,...(1)  \  .  f  \du]  (1)  , 


=  /  («• 
Jan  \ 


ds  - 


Lt  [dv\ 


=  Lm -  ds + L  (l^Ki  - “[fr1)  *  (3,7) 

=  L  *({*  ~  w’‘-' )  * + ^  L ^  (6  - 

To  obtain  the  last  identity  we  have  used  the  fact  that  (k  =  0  on  cr*.  and  that 

/[«*,*—/  ft*—/  ®=£t*-o. 

y<rfc  cw  ./an  at/  Jan  ai/ 

Concerning  the  two  terms  in  the  last  expression  of  equation  (3.17),  it  is  reasonable  to  assume 
that  the  first  term  /an  xj)((k  —  u^j)  ds  will  be  the  larger  (at  least  for  moderate  size  cracks). 
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If  we  now  substitute  i}>  =  8p  —  8q  into  the  last  expression  in  (3.17)  and  disregard  the  second 
term  we  get 

m*  *  (&  -  *&)(P)  ~  (&  -  «£!,)«)•  (3.18) 

We  now  chose  xf)<Tl  =  8p0  —  8p{  so  that  Pq  maximizes  (£1  —  u>^t)(P)  and  P\  minimizes 
(£1 —  wiiliXQ)-  For  subsequent  2  <  k  <  n  we  chose  i^ak  =  8p0  —  8pk  (the  same  Pq  as  for  k  =  1) 
where  Pk  is  selected  so  that  the  expression  (£jt  —  io^i)(Fb)  —  (€k  —  w^,i ){Pk)  is  °f  maximal 
magnitude,  subject  to  the  constraint  that  the  points  P0,  Pi,  Pi, . . . ,  Pn  stay  well  separated. 
We  cannot  allow  any  of  the  Pj  to  coincide,  for  then  we  would  be  duplicating  boundary 
measurements,  potentially  leading  to  a  singular  Jacobian  as  described  in  the  remark  before. 
The  above  construction  makes  the  two-electrode  currents  appear  somewhat  like  the  currents 
used  for  the  uniqueness  result  proven  in  [4].  For  the  uniqueness  result  we  needed  n  +  1 
two-electrode  currents  with  n  +  2  distinct  electrode  locations  -  the  currents  prescribed  above 
only  number  n  (with  n  +  1  electrode  locations).  We  expect  that  this  deficiency  is  more  than 
compensated  by  the  fact  that  the  electrode  locations  change  as  the  iterations  proceed. 


4  Integral  Equation  Formulation  and  Discretization 


We  now  proceed  to  formulate  the  boundary  value  problem  (3.2)  as  an  integral  equation  on 
the  boundary  of  the  region  Q  \  £.  Let  T(i,y)  denote  the  fundamental  solution  for  the  two 
dimensional  Laplacian  given  by 

r(X,y)  =  2^logdX  _yl)’  x’2/GlR2- 

The  application  of  standard  potential  theory  arguments  (see  [7],  Sections  3.B-3.D)  shows 
that  u,  the  solution  to  (3.2),  for  x  €  ft  \  E,  can  be  represented  as 

d  . .  A  r  .  raul 


«(*)  =  Ja n'U^!fr~r(X' y )  “  r^’  dSy  ~  'p  Lk  r^’ y )  ^ 


(y)dsy.  (4.1) 


Here^r-  on  dPl  denotes  the  normal  outward  derivative  with  respect  to  the  y  variable.  On 
any  of  the  cracks  <7 *,  v  denotes  a  unit  normal  field  and  |^|  =  denotes  the  jump 

in  the  normal  flux  across  the  crack.  To  simplify  notation  we  shall  use  the  notation  <f>k  for 
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the  jump  [|^j  across  ov  The  formula  (4.1)  expresses  the  value  of  u  at  any  point  in  ft  \  E  in 
terms  of  the  Dirichlet  and  Neumann  data  for  u  on  5ft  and  the  jump  in  across  the  cracks 
E.  This  is  simply  Green’s  third  identity  applied  to  the  region  ft  \  E.  It  is  straightforward 
to  check  that  since  has  at  most  an  r-1/2  singularity  at  the  endpoints  of  any  crack  and 
since  the  singularity  of  T  at  x  =  y  is  only  logarithmic,  equation  (4.1)  also  holds  for  x  €  E. 
Because  u  is  a  constant  c/  on  <7/  (and  continuous  in  ft)  this  implies  that 

/  (u(y)ir-r(x,y)-r(x,y)ilf(y))dsy-'jl2  f  r(x,y)My)  dsy  =  ct  (4.2) 
Jd O  OVy  Jok 

for  x  €  <Ji .  For  x  €  5ft,  an  argument  similar  to  that  which  led  to  (4.1)  leads  to  the  equation 

~lu(x)+  f  u(y)-S-r(x,y)dsy  -  £  f  T(x,y)<f>k(y)  dsy  =  f  T(x,y)rp(y)  dsy.  (4.3) 
l  Jd n  OVy  rr  Jak  J9(i 

As  discussed  earlier,  if  the  constants  c/  are  treated  as  unknowns  they  can  be  determined 
uniquely  from 

f  (f>ids  =  0,  /  =  l,...n,  (4.4) 

J&1 

with  the  normalization  fanuds  =  0.  Combining  these  n  +  1  conditions  with  equations  (4.2) 
and  (4.3)  we  arrive  at  the  following  system  of  integral  equations 


1  f  d  t 

-~u(x)+  u(y)—-T(x,y)dsy  -  ^{x,y)<f>k(y)  dsy 

£  Jd  17  OVy  j  J<rk 

=  /  T (x,y)i/>(y)dsy,  x  6  5ft 

Jd  n 


f  d 

Lu(v)a7?(x'v)ds’  - 


n  [ 

r(*>y)My)dsy  -  q 

k= i  J°k 

I  r(x,y)ip(y)dsy,  x  E  <7t 
Jdn 


with  the  additional  constraints 


/  fads  =  0,  /  =  1, . . . n, 

JO\ 

f  uds  =  0. 

Jdn 


The  unknowns  here  are  the  value  of  u  on  5ft,  <f>k  on  cr*.  and  the  constants  c/,  l  —  1, . . .  ,n. 
Given  a  set  of  cracks  E  and  Neumann  data  one  can  solve  equations  (4.5),  (4.6)  and  (4.7) 
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to  obtain  these  quantities.  The  solution  to  the  boundary  value  problem  (3.2)  at  any  point 
in  Cl  can  then  be  obtained  from  the  representation  (4.1). 

One  useful  fact  to  note  is  the  following.  Let  uo  be  a  harmonic  function  on  Cl  with 
Neumann  data  4>-  The  same  reasoning  used  to  derive  the  integral  equations  for  u  shows  that 
u0  satisfies  the  boundary  integral  equation 

\  t  Q  r 

-  2U°(*)  +  Jd(l  y)  dsv  =  Jgil  r(z> y)H y)  dsv  (4.8) 

for  x  €  dCl.  A  unique  solution  again  requires  a  normalization  such  as  /an  u0  ds  =  0.  If  x  6  Cl 
then  uo(x)  can  be  represented 

u°(x)  =  J^(u0(y)—r(x,  y )  -  T(x,  y)ip(y))  dsy.  (4.9) 

Let  v  denote  the  difference  u  —  «o-  Combining  equation  (4.5)  with  (4.8)  and  combining  (4.6) 
with  equation  (4.9),  we  find  that  v  satisfies 

~\v(x)  + Janv(y)'i^T(x’y)dsy-it, fa  T(x>y)My)dsy  =  xedci  (4.io) 

y  k—i 

I  v(y)-S-r(x,  y)  dsy-^2  f  r(x,  y)<t>k{y)  dsv  ~ci  =  -«o(*)>  x  € 

JdVl  Wy  JOk 

where  <f>k  denotes  the  jump  in  the  normal  derivative  j^v  across  <rk.  Note  that  this  is  the 
same  as  the  jump  in  ^ju,  since  u0  is  smooth  in  Cl.  The  conditions 


J'  <}>ids  =  0,  /  =  l,...n, 

<T\ 

/  v  ds  =  0, 

Jan 


(4.11) 


are  also  still  enforced.  Equations  (4.10)  and  (4.11)  provide  a  means  of  directly  computing 
the  perturbation  v  =  u  —  tt0  caused  by  the  presence  of  the  cracks.  This  formulation  is  more 
advantageous  than  the  original  formulation  since  we  will  use  singular  boundary  data  0.  The 
function  v  =  u  —  uq,  however,  is  smooth  up  to  dfi  and  hence  avoids  any  of  the  problems 
associated  with  the  lack  of  regularity  of  t/>.  Moreover,  for  the  specific  two-electrode  Neumann 
data  (and  a  domain  in  the  form  of  a  ball,  as  considered  later)  we  have  a  closed  form  solution 
for  u0. 

Having  derived  the  boundary  integral  formulation  we  now  briefly  discuss  how  we  discretize 
it  by  means  of  the  so-called  Nystrom  method.  Suppose  that  the  boundary  of  the  region  0 
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is  parameterized  by  z(t)  =  (z\(t),Z2(t)),  0  <  t  <  1.  Let  each  crack  crk  be  parameterized  by 
z(t)  for  k  <  t  <  k  +  1.  In  terms  of  this  parameterization,  equations  (4.10)  and  (4.11)  can  be 
written  as 

“2  +(*)+JoK(*>t)<K*)*-5lJk  G(s,t)<f>(t)dt  =  0,  s€[0,l)  (4.12) 

r1  JL  rk+l 

/  K(s,t)<j>(t)dt  -  /  G(s,t)<f>(t)dt  -  ci  =  -u0(z(s)),  s  €[/,/+ 1) 

Jo  k= i  •'* 

and 

/■H-l 

l  <f>{t)\z\t)\dt  =  0,  /  =  1,. . .  ,n, 

/  <^(i)|^r/(<) |  dt  =  0 
Jo 

where  K(s,t)  =  g^-r(^(s),z(t))|z,(t)j  and  G(s,<)  =  r(z(s),  z(f))|2r'(<)|.  The  functions  u  and 
<l>k  have  also  been  replaced  by  the  single  function  <j)  defined  on  [0,n  +  1]  by  <f>(t)  =  v(z(t )) 
for  t  €  [0, 1),  (j>(t )  =  <f>k{z{t)  for  t  €  [&,  k  +  1). 

Let  tj  and  Wj,  j  =  1, . . . ,  m,  denote  the  nodes  and  weights  of  a  quadrature  rule  on  [0, 1], 
so  that 

/  f(t)dt  « 

Jo  j=i 

for  reasonably  smooth  functions  f(t).  Nystrom’s  method  for  solving  an  integral  equation  of 
the  form 

4>{s)  +  [  K(s,  t)<f>(t)  dt  =  f(s ) 

Jo 

consists  of  replacing  the  integral  by  a  quadrature  rule  to  obtain 

m 

<!>{*)  +  ]C  tf(Mitofe)  =  f(s)- 

j= i 

Letting  s  assume  the  values  <i, . . . ,  tm  we  obtain  the  m  x  m  linear  system 

m 

<f>i  T  'y  ^  Kij  <f>j  =  fi  i  =  1, . . . ,  m 

3= 1 

where  <fn  =  <£(<,),  =  K{ti,tj)u3j  and  /;  =  /(t*).  The  intention  is  that  <f>i  =  <^(f;)  ss 

A  complete  treatment  of  Nystrom’s  method  for  second  kind  Fredholm  equations  can  be  found 
in  [2]. 


17 


For  a  first  kind  integral  equation  with  a  smooth  kernel, 

Tk4>=  ('  K(s,t)<f>(t)dt  =  f(s), 

Jo 

direct  discretization  usually  leads  to  a  linear  system  which  is  very  poorly  conditioned.  This 
stems  from  the  fact  that  the  corresponding  inverse  operator  T^1  is  unbounded  on  whatever 
space  the  equation  is  posed,  typically  <7(0, 1)  or  L2( 0, 1).  The  singular  values  of  the  operator 
Tk  approach  zero  rapidly,  so  that  the  solution  <f>  is  extremely  sensitive  with  respect  to  /  or  to 
noise  on  the  right  hand  side  of  the  equation.  The  smoother  the  kernel,  the  faster  the  singular 
values  decay  and  the  poorer  is  the  conditioning  of  the  linear,  system.  The  solution  of  first 
kind  integral  equations  therefore  often  requires  some  kind  of  regularization.  As  discussed 
in  [5],  however,  if  the  kernel  K(s,t)  is  singular  enough  on  the  diagonal  s  =  t,  reasonable 
results  can  be  obtained  from  a  direct  discretization  of  the  equation,  without  regularization. 
In  the  case  of  equations  (4.12)  the  first  kind  portion  of  the  equations  (on  the  cracks)  have 
a  logarithmic  singularity  along  the  diagonal  s  =  t,  and  hence  regularization  should  not  be 
necessary.  We  have  applied  Nystrom’s  method  directly  to  the  equations.  The  linear  systems 
obtained  in  this  manner  have  good  conditioning  and  in  all  test  cases  in  which  we  have  a 
closed  form  solution,  this  method  has  produced  solutions  in  complete  agreement  with  the 
closed  form  solution. 

To  apply  Nystrom’s  method  to  the  equations  (4.12),  we  replace  the  integrals  over  the 
intervals  [/,  /  +  1],  /  =  0, . . .  ,n,  with  the  quadrature  rule  and  then  let  s  assume  the  values 
/  +  ti,  i  =  l,...,m,  /  =  0, ...,n.  This  yields  the  following  linear  system  in  the  variables 

^01  j  •  •  •  >  *Ao m i  -  •  •  t  »  C| ,  .  .  .  ,  ^ . 


i  to  n  m 

+  ^2  K{ii,  tj)vj<f>oj  ^  ^  G(U,  k  +  tj)uj<f>kj 

1  J=l  fc=l j= 1 


Y,K(l  +  tj)u>j<f>0j  -'52'5l,G{1  +  ti,  k  -f-  tj)uj<i>kj  -  Cl 

j-l  k=l j= 1 

i  =  l,...,m,  / 


0,  (4.13) 

l,...,m 

-u0(z(l  +  ti)), 
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and 


'52<t>ijUj\z'(l  +  tj)\  =  0,  /  =  l,...,n, 

i= l 

m 

'£<t>°j“>\z'(tj)\  =  o, 

j- 1 

where  the  intention  is  that  4>kj  «  <f>(k+tj)  and  c*  «  c*.  In  this  formulation  there  are  actually 
mn  +  m  +  w  4- 1  equations  for  the  mn  +  m  +  n  unknowns  <f>kj,  k  =  0, . . . ,  n,  j  =  1 , . . . ,  m  and 
Cfc,  fc  =  1, . . .  n.  However,  a  careful  analysis  shows  that  the  first  m  equations  have  a  linear 
dependency  (the  coefficients  sum  to  zero)  so  that  any  one  of  them,  e.g.,  the  first,  can  be 
dropped.  This  gives  a  linear  system  of  mn  +  m  +  n  equations  for  the  unknowns  <f>kj  and  c*. 

One  need  not  use  the  same  quadrature  rule  on  the  boundary  of  0  and  the  cracks. 
Since  the  solution  is  smooth  on  dto,  we  allocate  evenly  spaced  nodes  there  by  t,  =  i/m, 
i  =  0, . . .  ,m  —  1.  The  weights  are  simply  w,  =  1/m,  corresponding  to  the  trapezoidal  rule 
on  the  closed  curve  dfl.  On  each  crack  <f>  will  typically  have  r-1^2  singularities  at  the  end¬ 
points  and  so  a  quadrature  rule  is  chosen  which  places  more  nodes  near  these  singularities 
(but  not  actually  at  the  endpoints).  If  the  linear  crack  with  endpoint  at  (a,  6),  angle  6  and 
length  A  is  parameterized  as  (a  +  tXcosO,  b+  <Asin0),  0  <  t  <  1,  then  the  nodes  are  chosen 
as  (assuming  m  is  even) 

r,i  —  1/2.  .  ,  m 

ti  —  /(  ~  )>  *  =  1 » *  *  * » » 

m  i 

and  ti  =  1  —  for  i  =  y  +  1,...  ,m,  where  /( x)  =  29~1xq  and  q  is  a  positive  number. 

The  parameter  q  controls  the  spacing  of  the  nodes  with  q  >  1  causing  the  nodes  to  “bunch 
up”  near  0  and  I.  We  have  used  q  =  2.5.  The  weights  w,  are  chosen  as 


HU  +  *2)7 


i  =  l 


Ui  ~  '  §(<.+1  “  t  =  2, . . . , m  —  1 

1  2(^*7*— 1  T  ) j  t  =  m 


corresponding  to  a  midpoint  rule  with  variably  spaced  nodes.  We  have  used  60  nodes  on 
both  dQ  and  each  crack. 

One  difficulty  with  discretizing  equations  (4.12)  is  the  presence  of  the  logarithmic  singu¬ 
larity  in  the  first  kind  portion  of  the  equations  along  the  diagonal  s  =  t.  We  deal  with  this 
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by  using  a  simple  form  of  product  integration  based  on  our  quadrature  rule.  See  [2],  Section 
3.2,  for  more  details. 


5  The  Jacobian  and  Newton’s  Method 

We  recall  that  a  central  part  of  our  crack  reconstruction  algorithm  is  to  find  a  solution  of 
the  4 n  x  4n  system 

G*(E)  =  F(E,^,wffJ-F (5.1) 

i 

=  F(E,0„fc,w<rJ  —  f  A:  =  1, ... ,n. 


Here  the  four  components  of  F(E,V>,  w)  are  given  by 

...  r 

<’>)=/  u(E,V)^-<fo  i=l . 4, 

Jan  av 

where  u  solves  (3.2)  and  are  the  functions  from  (3.9)-(3.11).  If  we  use  the  vector  Q  6  IR4n 
to  describe  the  crack  configuration  E  and  we  parameterize  9(0  \  E)  as  described  in  the  last 
section,  then  the  discretized  version  of  each  of  these  components  becomes 

m  a 

HQ^iw)  =  Y,uj-fo{ti)\z\ti)\ui  (5.2) 

where  tj  is  the  jth  node  for  Nystom’s  method  on  90,  uij  is  the  corresponding  weight.  The 
variables  u}  are  given  by  u,  =  uo(tj)  +  <po j,  where  <f>oj  form  the  first  part  of  the  solution  to 
the  system  (4.13). 

Our  approximate  method  for  the  solution  of  (5.1)  is  a  variation  of  Newton’s  method.  For 
that  we  need  an  effective  method  for  the  calculation  of  the  Jacobian.  If  q  denotes  one  of  the 
components  of  Q ,  differentiation  of  the  expression  (5.2)  yields 

9F  ^  ( duj  dw/A  9  (dw,s^\\  oX 

Here  we  have  assumed  that  \z'\  is  independent  of  q  on  9f l  and,  as  mentioned  previously, 
we  have  ignored  the  functional  dependence  of  the  applied  current  flux  ip  on  E  (this  in 
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particular  gives  that  =  0).  Note  that  the  functions  w  defined  by  (3.9)-  (3.11)  are  indeed 
differentiable  with  respect  to  the  parameters  describing  the  cracks. 

To  evaluate  ^  from  equation  (5.3)  we  therefore  need  to  calculate  |^,  the  derivatives  of 
the  solution  to  (4.13)  with  respect  to  the  parameters  which  describe  the  cracks.  These  can 
be  computed  with  little  additional  effort.  Let  the  linear  system  (4.13)  be  written  in  the  form 

A(Q)<j>(Q)  =  f(Q)  (5.4) 

where  A(Q)  is  the  mn  +  m  +  n  by  mn  +  m  +  n  matrix  appearing  in  (4.13),  f{Q)  €  iRmn+m+" 
is  the  right  hand  side,  and  <f>(Q )  6  JR,mn+Tn+n  }s  the  solution  to  the  system  (including  the 
constants  c).  Of  course  all  of  these  quantities  depend  on  the  parameters  Q.  Differentiation 
of  equation  (5.4)  with  respect  to  any  one  of  the  parameters  of  Q  and  use  of  the  fact  that 

=  0  gives 

A^frfq-aiqm=Jim'  (5-5) 

so  the  derivative  satisfies  a  linear  equation  of  exactly  the  same  form  as  <f>  but  with  a  different 
right  hand  side.  Once  (5.4)  has  been  solved  for  <j>(Q)  (e.g.,  by  LU  decomposition),  equation 
(5.5)  can  be  solved  for  by  simply  computing  the  right  hand  side  and  reusing  the  LU 
decomposition. 

Below  is  a  global  description  of  our  algorithm  for  crack  reconstruction.  We  denote  by 
E*  =  { cr* } x  the  estimated  cracks  at  the  ith  stage  of  the  algorithm,  and  by  Q'  we  denote 
the  corresponding  set  of  parameters. 

1.  Make  an  initial  guess  E°,  set  i  —  0. 

2.  Select  the  maximally  sensitive  two-electrode  fluxes  Vvj.  corresponding  to  the  cracks  u'k, 
k  =  1, . . .  ,n,  in  the  sense  defined  in  section  3. 

3.  Measure  (simulate)  the  boundary  voltage  data  for  each  flux  and  compute  f(^<Tj_,  w^), 
*  =  l,...,n. 

4.  Compute  the  voltage  data  for  each  flux  at  E  =  E’  and  compute  F(E’, 
k  =  l,...,n. 
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5.  Compute  Gjt(S’)  =  F(E‘,  w^)  —  wCTi ),  k  =  1  Let  G(E)  denote  the 

vector  of  all  residuals{Gjfc(E)}£=1  (a  4n-vector).  If  G(E‘)  (=  G(Q'))  is  sufficiently 
small,  terminate  with  the  answer  E  =  E*. 

6.  Compute  the  approximate  Jacobian  J,  =  “Z)qG((5)|q=q«”,  in  the  sense  described  ear¬ 
lier. 

7.  Compute  the  Newton  update  8Q'  by  solving  J{8Q'  —  —G(QX). 

8.  Update  Qt+1  =  Q'  +  8Q ',  *  =  i  +  1,  and  go  to  step  2. 

In  our  implementation  of  Newton’s  method,  as  in  [15],  we  solve  the  linear  system 

J,8Q'  =  -G(Q‘) 

subject  to  the  constraint 

II  W)ll  <  P 

where  A  is  an  appropriately  chosen  diagonal  weighting  matrix  and  p  a  specified  parame¬ 
ter.  This  constrained  problem  is  solved  in  the  least  squares  sense  by  means  of  a  Lagrange 
multiplier  method  as  outlined  in  [11].  The  constraint  on  the  update  markedly  improves  the 
behavior  of  Newton’s  method  far  from  the  solution. 

It  should  be  mentioned  that  we  also  constrain  the  cracks  to  lie  inside  the  domain  U; 
if  the  algorithm  attempts  to  move  a  crack  outside  f!  we  perform  a  simple  reflection  of  the 
endpoint(s)  that  lie  outside  to  get  get  the  updated  crack  fully  inside  Q.  The  algorithm  very 
rarely  attempts  to  move  a  crack  outside  the  domain  (unless  the  data  based  on  which  we 
perform  the  reconstruction  is  generated  from  a  collection  of  cracks,  at  least  which  of  one  lies 
very  near  the  boundary).  Our  implementation  will  also  allow  the  cracks  to  intersect  each 
other;  the  integral  equation  formulation  continues  to  make  sense  in  this  case,  although  in 
practice  one  must  be  careful  when  dealing  with  the  logarithmic  singularities  which  arise  as 
a  result  of  the  intersection.  We  deal  with  this  by  again  using  a  form  of  product  integration. 
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6  Computational  Experiments 


We  have  performed  a  significant  amount  of  computational  experimentation  with  the  algo¬ 
rithm  we  have  developed.  In  this  section  we  briefly  describe  five  such  experiments  which 
we  find  to  be  representative.  The  experiments  are  all  performed  on  simulated  data.  It  is 
our  goal  to  eventually  try  our  algorithm  on  actual  experimental  data.  A  project  to  build 
an  experimental  setup  and  a  data  gathering  device  is  currently  in  progress  [10].  To  go  part 
of  the  way  towards  reconstruction  based  on  real  data,  the  last  of  the  experiments  described 
here  pertains  to  simulated  data  with  a  built-in  level  of  noise. 

The  domain  on  which  the  reconstruction  is  performed  is  in  all  cases  the  unit  disk.  The 
graphics  by  means  of  which  we  illustrate  the  progress  of  the  algorithm  is  the  same  for  all 
experiments.  Each  step  in  the  iteration  is  illustrated  by  a  picture  containing  two  copies  of 
the  unit  disk.  The  disk  on  the  left  depicts  the  previously  estimated  locations  of  the  cracks 
(a  set  of  line  segments)  as  well  as  the  boundary  locations  of  the  optimally  sensitive  electrode 
locations.  The  electrode  locations  are  marked  with  small  circles;  the  circle  corresponding  to 
Pq  has  been  darkened.  It  is  voltage  data  corresponding  to  the  currents  generated  by  these 
electrode  locations  that  are  used  for  the  iterative  update.  The  updated  estimates  of  the 
crack  locations  are  shown  as  solid  line  segments  in  the  disk  on  the  right.  The  true  cracks 
which  the  algorithm  seeks  to  reconstruct,  i.e.,  the  cracks  from  which  the  simulated  boundary 
voltage  data  is  generated,  are  shown  as  dashed  line  segments  (or  dashed  circular  arcs)  in  the 
disk  on  the  right. 

In  our  first  experiment  the  simulated  data  comes  from  three  cracks,  two  of  which  are 
very  near  the  boundary.  They  are  each  0.05  units  away  from  the  boundary.  The  cracks  have 
lengths  0.25,  0.3  and  0.35.  We  start  the  reconstruction  procedure  by  attempting  to  fit  the 
simulated  data  with  data  generated  from  a  single  crack.  As  the  initial  guess  we  select  a  crack 
joining  the  points  (-0.2,  0.2)  to  (-0.2,  0.6).  Figure  la  shows  the  first  step  of  the  iteiation. 
After  31  steps  the  algorithm  has  found  a  root  for  the  system  of  four  functionals  and  reduced 
the  residual  error  (|G(Q)|)  to  8.94  x  10-13.  Figure  lb  shows  step  31  of  the  iteration. 
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Figure  lg:  Iteration  19:  four  cracks  fitted  to  three  crack  boundary  data. 


At  this  point  it  is  not  possible,  based  on  four  functionals,  to  determine  whether  the 
simulated  boundary  voltage  data  comes  from  just  a  single  crack  or  several  more.  The  only 
way  to  determine  this  is  to  try  to  fit  more  functionals.  We  take  the  crack  shown  in  the  right 
disk  in  Figure  lb  and  divide  it  into  two  cracks  by  cutting  out  a  piece  1/10  of  the  length  at 
the  center.  The  two  resulting  cracks  are  now  used  as  the  initial  guess  for  our  algorithm  based 
on  8  functionals.  The  first  step  of  the  two-crack  iteration  is  show  in  Figure  lc.  After  18 
steps  the  “two-crack  iteration”  has  found  a  root  to  the  8- variable  system  and  has  reduced  the 
residual  error  to  6.9  x  10-16.  The  final  step  is  shown  in  Figure  Id.  Finally  we  take  the  largest 
of  the  two  cracks  from  the  right  disk  in  Figure  Id,  divide  it  into  two  pieces  (by  cutting  out 
the  middle  1/10),  and  give  the  resulting  three  cracks  as  initial  guess  to  our  algorithm  based 
on  12  functionals.  This  “three  crack  iteration”  locates  the  root  after  24  steps,  reducing  the 
residual  error  to  2.0  x  10“n.  Steps  10  and  24  of  this  iteration  are  shown  in  Figure  le  and 
Figure  If,  respectively.  If  at  this  point  we  take  and  divide  one  of  the  three  cracks  again  and 
give  the  resulting  four  cracks  as  an  initial  guess  to  our  algorithm  based  on  16  functionals, 
then  one  of  two  things  is  likely  to  happen:  1)  these  new  four  crack  will  remain  essentially 
where  the  three  cracks  already  are  (and  the  residual  will  be  very  small)  2)  the  algorithm  will 
shrink  one  of  the  cracks  to  zero  length  and  the  three  remaining  will  stay  as  before.  In  the 
latter  case  the  residual  will  not  become  quite  as  small  since  we  impose  a  lower  limit  (of  0.01) 
on  the  length  of  the  cracks.  In  both  cases  the  behavior  clearly  indicates  that  three  cracks 
are  the  right  number  to  fit  the  data  (also  for  16  functionals).  Figure  lg  shows  step  19  in  a 
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“four  crack  iteration”  where  the  largest  crack  in  Figure  If  has  been  divided  into  two  pieces. 
As  is  evident  the  latter  of  the  two  possibilities  from  before  emerges  (the  residual  after  19 
steps  is  5.2  x  10-4). 

The  strategy  for  gradually  increasing  the  number  of  cracks  as  outlined  above  has  emerged 
after  a  significant  amount  of  numerical  experimentation  with  many  sets  of  simulated  data. 
The  alternate  strategy:  to  initially  guess  a  sufficient  number  of  cracks  and  then  let  the  itera¬ 
tions  proceed  to  convergence  has  generally  shown  itself  to  lead  to  much  slower  convergence. 
Our  method  based  on  the  use  of  only  a  few  functionals  and  the  adaptive  movement  of  elec¬ 
trodes  has  in  many  of  our  experiments  proven  itself  to  be  superior  to  a  least  squares  fitting 
algorithm  (using  the  entire  set  of  boundary  voltage  data,  but  a  fixed  set  of  electrodes).  Our 
experience  with  this  least  squares  approach  has  been  that,  except  for  the  one  crack  case, 
it  requires  a  very  accurate  initial  guess  in  order  to  converge  to  the  correct  solution.  For 
multiple  cracks  the  least  squares  functional  appears  to  contain  many  local  minima. 

This  is  not  to  say  that  our  approach  may  not  occasionally  be  somewhat  slower  than  indi¬ 
cated  by  the  previous  example.  We  illustrate  this  with  a  reconstruction  based  on  simulated 
data  from  four  cracks.  Figure  2a  shows  the  first  step  using  a  single  crack  (four  function¬ 
als).  After  17  steps  our  algorithm  finds  a  root  with  the  residual  reduced  to  1.72  x  10~14. 
However,  as  seen  in  Figure  2b  (which  shows  the  final  iteration)  the  single  crack  that  is  con¬ 
sistent  with  the  four  functionals  does  not  lie  near  any  of  the  four  cracks  that  were  used 
to  generate  the  data.  We  now  divide  this  single  crack  into  two  by  cutting  out  the  middle 
1/10.  Using  the  resulting  two  cracks  as  initial  guess  our  algorithm  based  on  8  functionals 
now  takes  a  considerable  number  of  steps  before  the  residual  is  even  reasonably  small  (and 
the  cracks  are  of  reasonable  size).  Figures  2c  and  2d  show  iterations  99  and  198,  respec¬ 
tively.  We  take  the  two  cracks  after  198  iterations  and  divide  each  of  them  into  two  using 
the  same  method  as  before.  The  resulting  four  cracks  are  provided  as  initial  guess  for  our 
algorithm  based  on  16  functionals.  Iterations  50  and  115  of  this  process  are  shown  in  Fig¬ 
ures  2e  and  2f,  respectively.  After  127  iterations  a  root  is  found  and  the  residual  has  been 
reduced  to  1.65  x  10~14.  Even  though  the  algorithm  is  extremely  slow  it  does  ultimately 
converge.  It  would  require  an  extremely  accurate  initial  guess  to  get  the  least  squares  algo¬ 
rithm  we  described  before  to  converge.  For  comparison  figure  2g  shows  the  eleventh  and  final 
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Figure  2g:  Least  squares,  iteration  11:  four  cracks  fitted  to  four  crack  data. 


iterate  of  a  least  squares  approach.  This  computation  was  done  using  the  full  boundary 
data  and  a  Levenberg-Marquardt  algorithm  as  outline  in  [11].  The  initial  guess  used  (con¬ 
sisting  of  four  cracks)  is  the  same  as  that  used  for  figures  2e  and  2f.  The  least  squares 
approach  quickly  locates  a  local  minimum  and  terminates.  The  four  cracks  that  are  used 
appear  to  merge  into  two. 

Frequently  in  practice  there  are  either  a  fairly  limited  number  of  well-separated  cracks, 
or  many  cracks  clustered  in  certain  locations.  Figures  3a  and  3b  show  the  first  iteration 
and  the  fifth  iteration  using  our  algorithm  with  one  crack  (and  four  functionals)  to  fit  simu¬ 
lated  data  coming  from  10  cracks  located  in  two  clusters.  A  root  is  found  at  the  fifth  iteration 


Figure  3a:  Iteration  1 :  one  crack  fitted  to  ten  crack  boundary  data. 
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with  residual  7.11  x  10-15.  We  now  cut  this  last  crack  in  two  (deleting  the  middle  1/10)  and 
give  the  resulting  two  cracks  as  input  to  our  algorithm  based  on  8  functionals.  The  algorithm 
locates  a  root  after  15  iterations,  the  corresponding  crack  locations  are  shown  in  figure  3c. 
We  again  split  the  larger  of  the  two  cracks  and  use  that  as  input  for  a  “three  crack”  version 
of  our  algorithm  -  the  result  after  99  steps  is  shown  in  figure  3d,  although  the  algorithm  has 
not  yet  located  a  root. 

Another  situation  that  would  occur  in  practice  is  when  the  “real”  cracks  have  shapes 
other  than  line  segments.  Figure  4a  and  4b  show  the  first  and  twelfth  iteration  of  our  algo¬ 
rithm  using  one  crack  in  a  case  where  the  simulated  data  is  generated  by  a  crack  in  the  shape 


Figure  4a:  Iteration  1:  one  crack  fitted  to  curved  crack  boundary  data. 


Figure  4b:  Iteration  12:  one  crack  fitted  to  curved  crack  boundary  data. 
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of  a  circular  arc  centered  at  the  origin  and  going  from  (0.0,0. 5)  to  (-0.5, 0.0).  The  circu¬ 
lar  arc  is  indicated  by  the  dashed  lines  in  the  disks  on  the  right  -  these  are  not  individual 
cracks  but  are  used  to  indicate  the  curve.  The  algorithm  finds  a  root  at  the  twelfth  iteration. 
We  now  proceed  to  divide  the  last  crack  of  the  “one  crack  iteration”  into  two  (as  before)  and 
we  iterate  using  the  “two  crack  version”  of  our  algorithm.  Figure  4c  shows  the  50th  iteration 
of  this  process.  The  algorithm  has  not  yet  located  a  root  -  actually  the  crack  locations  have 
not  changed  considerably  since  the  15th  iteration. 


Figure  4d:  Iteration  50:  four  cracks  fitted  to  curved  crack  boundary  data. 


We  now  divide  the  final  two  cracks  from  Figure  4c  into  four  and  give  these  as  an  initial 
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guess  to  the  “four  crack  version”  of  our  algorithm.  Figure  4d  shows  iteration  50  using  four 
cracks.  The  algorithm  has  not  yet  located  a  root.  The  two  longest  cracks,  which  already 
emerge  after  13  iterations,  provide  a  reasonable  approximation  to  the  curved  crack.  One  of 
the  two  shorter  cracks  stays  close  to  the  circular  arc  (and  the  two  large  cracks)  the  other 
one  seeks  to  exit  the  domain;  the  program  apparently  cannot  adjust  them  to  find  a  root,  but 
there  is  a  lower  bound  on  their  length,  so  they  cannot  entirely  disappear. 

In  the  final  example  we  have  taken  data  generated  by  a  two  cracks  and  added  10% 
noise.  The  noise  added  is  independent  and  gaussian  with  a  zero  mean  and  standard  de¬ 
viation  equal  to  1/10  the  mean  square  value  of  u  —  uo  on  dCl,  where  u  is  the  potential 


Figure  5a:  Iteration  1:  two  crack  fitted  to  two  crack  data,  10  percent  noise. 


Figure  5b:  Iteration  24:  two  crack  fitted  to  two  crack  data,  10  percent  noise. 
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and  uo  is  the  harmonic  function  with  the  same  flux  as  u.  In  figure  5a  we  show  the  first 
iteration.  Figure  5b  shows  the  24th  iteration,  at  which  point  the  algorithm  terminates  hav¬ 
ing  found  a  root.  Figure  5c  shows  the  result  of  taking  the  final  crack  locations  in  figure 
5b  and  using  them  as  an  initial  guess  to  a  least-squares  minimization  routine.  This  rou¬ 
tine  uses  a  fixed  set  of  electrodes,  those  from  the  last  iteration  of  our  algorithm.  These 
are  presumably  close  to  the  most  sensitive  electrode  locations  for  the  given  cracks.  The 
Levenberg-Marquardt  routine  reduces  the  total  least-squares  residual  from  0.085  to  0.032  in 
13  iterations  and  improves  the  estimate  of  the  crack  locations. 


Figure  5c:  Least  squares,  iteration  13:  two  crack  fitted  to  two  crack  data,  10 
percent  noise. 


7  Summary 

In  this  paper  we  have  developed  a  very  efficient  algorithm  for  the  reconstruction  of  a  collec¬ 
tion  of  cracks  based  on  electrostatic  boundary  measurements.  We  use  a  “dual”  variational 
formulation  for  the  forward  electrostatic  problem,  and  solve  this  numerically  by  means  of  a 
Nystrom’s  approximation  of  the  corresponding  boundary  integral  equations.  Our  reconstruc¬ 
tion  is  based  on  adaptively  changing  the  current  patterns,  so  as  to  maximize  the  sensitivity 
of  the  measured  voltage  differences.  Our  reconstruction  is  based  on  a  limited  set  of  aver¬ 
ages  of  the  boundary  voltage  measurements  as  opposed  to  the  entire  set  of  measurements; 
this  should  lead  to  greater  efficiency  and  less  rigidity.  The  algorithm  is  currently  entirely 
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two-dimensional,  but  it  should  be  very  interesting  to  extend  it  to  three  dimensions.  At  this 
point  we  have  only  investigated  the  behavior  of  the  algorithm  when  used  on  “synthetic” 
data  (including  data  with  noise).  It  should  be  very  interesting  to  apply  the  algorithm  to 
data  coming  from  “real”  experiments.  Frequently  cracks  appear  as  clusters  of  many  small 
(microscopic)  cracks;  our  algorithm,  when  applied  to  data  generated  by  clusters  of  small 
cracks,  often  very  successfully  locates  a  set  consisting  of  a  few,  well-separated  cracks.  In  this 
context  it  should  be  extremely  interesting  to  analyze  in  what  sense  this  reflects  the  behavior 
of  the  forward  problem.  To  be  more  specific:  it  should  be  interesting  to  study  in  what  sense 
a  cluster  of  small  (microscopic)  cracks  in  an  appropriate  limit  approaches  a  single  (lumped) 
macroscopic  crack. 


36 


References 


[1]  Alessandrini,  G.,  Stable  determination  of  a  crack  from  boundary  measurements.  To 
appear,  Proc.  Roy.  Soc.  Edinburgh. 

[2]  Atkinson,  K.E.,  A  survey  of  numerical  methods  for  the  solution  of  fredholm  integral 
equations  of  the  second  kind ,  SIAM,  Philadelphia,  PA,  1976. 

[3]  Barber,  D.  and  Brown,  B.,  Recent  developments  in  applied  potential  tomography  -  APT, 
in  Information  Processing  in  Medical  Imaging,  Bacharach,  S.  ed.,  Nijhoff,  Amsterdam, 
1986,  pp.  106-121. 

[4]  Bryan,  K.  and  Vogelius,  M.,  A  uniqueness  result  concerning  the  identification  of  a 
collection  of  cracks  from  finitely  many  electrostatic  boundary  measurements.  To  appear, 
SIAM  J.  Math.  Anal. 

[5]  de  Hoog,  F.R.,  Review  of  fredholm  equations  of  the  first  kind.  In  the  application  and 
numerical  solution  of  integral  equations ,  Andersen,  R.S.  et  al  (eds.),  Sijthoff  and  No- 
ordhoff:  Alphen  aan  den  Rijn,  1980. 

[6]  Eggleston,  M.R.,  Schwabe,  R.J.,  Isaacson,  D.,  Goble,  J.C.  and  Coffin,  L.F.,  Three- 
dimensional  defect  imaging  with  electric  current  computed  tomography.  GE  Technical 
Report  91CRD039,  Schenectady,  NY. 

[7]  Folland,  G.B.,  Introduction  to  partial  differential  equations.  Princeton,  NJ:  Princeton 
University  Press,  1976. 

[8]  Friedman,  A.  and  Vogelius,  M.,  Determining  cracks  by  boundary  measurements.  Indiana 
Univ.  Math.  J.,  38  (1989),  pp  527-556. 

[9]  Gisser,  D.G.,  Isaacson,  D.  and  Newell,  J.C.,  Electric  current  computed  tomography  and 
eigenvalues.  SIAM  J.  Appl.  Math.,  50  (1990),  pp  1623-1634. 

[10]  Liepa,  V.,  Santosa,  F.  and  Vogelius,  M.  Crack  determination  from  boundary  mea¬ 
surements:  Computational  reconstruction  from  laboratory  measurements.  Submitted 
to  Journal  of  Nondestructive  Evaluation. 


37 


[11]  More,  J.,  The  Levenberg-Marquardt  algorithm:  implementation  and  theory.  Numerical 
Analysis  (Edited  by  Watson,  G.A.),  pp.  105-116.  Lecture  Notes  in  Math.  630.  Springer 
Verlag,  1977. 

[12]  Nedelec,  J.C.,  Integral  equations  with  non  integrable  kernels.  Integral  Eq.  Operator  Th., 
5  (1982),  pp  562-572. 

[13]  Nishimura,  N.,  Regularized  integral  equations  for  crack  shape  determination  problems. 
Inverse  Problems  in  Engineering  Sciences  (Edited  by  Yamaguti,  M.  et  al),  pp.  59-65. 
Springer  Verlag,  1991. 

[14]  Nishimura,  N.  and  Kobayashi,  S.,  A  boundary  integral  equation  method  for  an  inverse 
problem  related  to  crack  detection.  Int.  J.  Num.  Meth.  Eng.,  32  (1991),  pp  1371-1387. 

[15]  Santosa,  F.  and  Vogelius,  M.,  A  computational  algorithm  to  determine  cracks  from 
electrostatic  boundary  measurements.  Int.  J.  Eng.  Sci.  29  (1991),  pp  917-937. 

[16]  Yorkey,  T.,  Webster,  J.  and  Tompkins,  W.,  Comparing  reconstruction  algorithms  for 
electrical  impedance  tomography.  IEEE  Trans.  Biomedical  Eng.,  BME-34  (1987),  pp. 
843-852. 


38 


REPORT  DOCUMENTATION  PAGE 

Form  Approved 

OMB  No.  0704-0188 

PuWk:  reporting  burden  for  thi*  collection  of  information  «s  estimated  to  average  i  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this 
collection  of  information,  including  suggestions  for  reducing  this  burden  to  Washington  Headquarters  Services,  Directorate  for  information  Operations  and  Reports.  1215  Jefferson 
Oevrt  Highway.  Suite  1204.  Arlington,  V A  22202-4302.  and  to  the  Office  of  Management  and  Budget.  Paperwork  Reduction  Project  (0704-0 188).  Washington,  DC  20503 

1.  AGENCY  use  ONLY  (Leave  blank)  2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

June  1992  Contractor  Report 

4.  TITLE  ANO  SUBTITLE 

A  COMPUTATIONAL  ALGORITHM  FOR  CRACK  DETERMINATION. 

THE  MULTIPLE  CRACK  CASE 

S.  FUNDING  NUMBERS 

C  NAS1-18605 

C  NAS1-19480 

WU  505-90-52-01 

t.  AUTHOR(S) 

Kurt  Bryan 

Michael  Vogellus 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADORESS(ES) 

Institute  for  Computer  Applications  In  Science 
and  Engineering 

Mail  Stop  132C,  NASA  Langley  Research  Center 

Hampton,  VA  23665-5225 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

ICASE  Report  No.  92-24 

9.  SPONSORING /MONITORING  AGENCY  NAME(S)  ANO  AOORESS(ES) 

National  Aeronautics  and  Space  Adminstration 

Langley  Research  Center 

Hampton,  VA  23665-5225 

10.  SPONSORING /MONITORING 

AGENCY  REPORT  NUMBER 

NASA  CR-189665 

ICASE  Report  No.  92-24 

11.  supplementary  notes  Submitted  to  International  Journal 

Langley  Technical  Monitor:  Michael  F.  Card  of  Engineering  Sciences 

Final  Report 

1..  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

Uru-  <  ossified  -  Unlimited 

Subject  Category  64 

12b.  DISTRIBUTION  CODE 

13.  ABSTRACT  (Maximum  200  words) 

This  paper  develops  an  algorithm  for  recovering  a  collection  a  linear  cracks  in  a 
homogeneous  electrical  conductor  from  boundary  measurements  of  voltages  induced  by 
specified  current  fluxes.  The  technique  is  a  variation  of  Newton's  method  and  is 
based  on  taking  weighted  averages  of  the  boundary  data.  The  method  also  adaptively 
changes  the  applied  current  flux  at  each  iteration  to  maintain  maximum  sensitivity 
to  the  estimated  locations  of  the  cracks. 

14.  SUBJECT  TERMS 

inverse  problems;  nondestructive  evaluation;  boundary  integral 
equations 

Uaiig 

■ 

IS.  SECURITY  CLASSIFICATION 

OF  THIS  PAGE 

Unclassified 

19.  SECURITY  CLASSIFICATION 
OF  ABSTRACT 

20.  LIMITATION  OF  ABSTRACT 

N5N  7540*01-280-5500  Standard  Form  298  (Rev  2-89) 


Pretended  by  ansi  Std  239-18 
298-102 


NAS  A -Langley,  1992 


